Contracting cube#
In this demo we simulate a unit cube that is fixed at \(x = 0\) and free at \(x = 1\). We use a transversally isotropic material with fiber oriented in the \(x\)-direction.
import dolfin
try:
from dolfin_adjoint import (
Constant,
DirichletBC,
Expression,
UnitCubeMesh,
interpolate,
Mesh,
)
except ImportError:
from dolfin import (
Constant,
DirichletBC,
interpolate,
Expression,
UnitCubeMesh,
Mesh,
)
import pulse
from fenics_plotly import plot
pulse.iterate.logger.setLevel(10)
# Create mesh
N = 6
mesh = UnitCubeMesh(N, N, N)
# Create subdomains
class Free(dolfin.SubDomain):
def inside(self, x, on_boundary):
return x[0] > (1.0 - dolfin.DOLFIN_EPS) and on_boundary
class Fixed(dolfin.SubDomain):
def inside(self, x, on_boundary):
return x[0] < dolfin.DOLFIN_EPS and on_boundary
# Create a facet fuction in order to mark the subdomains
ffun = dolfin.MeshFunction("size_t", mesh, 2)
ffun.set_all(0)
# Mark the first subdomain with value 1
fixed = Fixed()
fixed_marker = 1
fixed.mark(ffun, fixed_marker)
# Mark the second subdomain with value 2
free = Free()
free_marker = 2
free.mark(ffun, free_marker)
# Create a cell function (but we are not using it)
cfun = dolfin.MeshFunction("size_t", mesh, 3)
cfun.set_all(0)
# Collect the functions containing the markers
marker_functions = pulse.MarkerFunctions(ffun=ffun, cfun=cfun)
# Create mictrotructure
V_f = dolfin.VectorFunctionSpace(mesh, "CG", 1)
# Fibers
f0 = interpolate(Expression(("1.0", "0.0", "0.0"), degree=1), V_f)
# Sheets
s0 = interpolate(Expression(("0.0", "1.0", "0.0"), degree=1), V_f)
# Fiber-sheet normal
n0 = interpolate(Expression(("0.0", "0.0", "1.0"), degree=1), V_f)
# Collect the mictrotructure
microstructure = pulse.Microstructure(f0=f0, s0=s0, n0=n0)
# Create the geometry
geometry = pulse.Geometry(
mesh=mesh,
marker_functions=marker_functions,
microstructure=microstructure,
)
# Use the default material parameters
material_parameters = pulse.HolzapfelOgden.default_parameters()
# Select model for active contraction
active_model = pulse.ActiveModels.active_strain
# active_model = "active_stress"
# Set the activation
activation = Constant(0.0)
# Create material
material = pulse.HolzapfelOgden(
active_model=active_model,
parameters=material_parameters,
activation=activation,
)
# Make Dirichlet boundary conditions
def dirichlet_bc(W):
V = W if W.sub(0).num_sub_spaces() == 0 else W.sub(0)
return DirichletBC(V, Constant((0.0, 0.0, 0.0)), fixed)
# Make Neumann boundary conditions
neumann_bc = pulse.NeumannBC(traction=Constant(0.0), marker=free_marker)
# Collect Boundary Conditions
bcs = pulse.BoundaryConditions(dirichlet=(dirichlet_bc,), neumann=(neumann_bc,))
# Create problem
problem = pulse.MechanicsProblem(geometry, material, bcs)
# Solve problem
pulse.iterate.iterate(problem, activation, 0.1)
2023-10-28 16:38:36,245 [1229] DEBUG pulse.iterate: Control: [0.0]
2023-10-28 16:38:36,246 [1229] DEBUG pulse.iterate: Target: [0.1]
2023-10-28 16:38:36,248 [1229] DEBUG pulse.iterate: Intial number of steps: 5 with step size 0.02
2023-10-28 16:38:36,248 [1229] INFO pulse.iterate: Iterating to target control (f_21)...
2023-10-28 16:38:36,249 [1229] INFO pulse.iterate: Current control: f_21 = 0.000
2023-10-28 16:38:36,250 [1229] INFO pulse.iterate: Target: 0.100
2023-10-28 16:38:36,251 [1229] DEBUG pulse.iterate: Check target reached: NO
2023-10-28 16:38:36,254 [1229] DEBUG pulse.iterate: Maximum difference: 1.000e-01
2023-10-28 16:38:36,255 [1229] DEBUG pulse.iterate: Try new control
2023-10-28 16:38:36,256 [1229] DEBUG pulse.iterate: Current control: 0.020
2023-10-28 16:38:36,430 [1229] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2023-10-28 16:38:36,431 [1229] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2023-10-28 16:38:36,432 [1229] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2023-10-28 16:38:36,433 [1229] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2023-10-28 16:38:36,433 [1229] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2023-10-28 16:38:36,434 [1229] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2023-10-28 16:38:36,435 [1229] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2023-10-28 16:38:36,436 [1229] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2023-10-28 16:38:36,437 [1229] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2023-10-28 16:38:36,438 [1229] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2023-10-28 16:38:36,438 [1229] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2023-10-28 16:38:37,926 [1229] DEBUG UFL_LEGACY: Blocks of each mode:
99 full
2023-10-28 16:38:38,257 [1229] DEBUG UFL_LEGACY: Blocks of each mode:
18 full
2023-10-28 16:40:26,172 [1229] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2023-10-28 16:40:26,173 [1229] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2023-10-28 16:40:26,174 [1229] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2023-10-28 16:40:26,175 [1229] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2023-10-28 16:40:26,176 [1229] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2023-10-28 16:40:26,177 [1229] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2023-10-28 16:40:26,178 [1229] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2023-10-28 16:40:26,178 [1229] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2023-10-28 16:40:26,179 [1229] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2023-10-28 16:40:26,180 [1229] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2023-10-28 16:40:26,181 [1229] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2023-10-28 16:40:26,439 [1229] DEBUG UFL_LEGACY: Blocks of each mode:
10 full
2023-10-28 16:40:26,738 [1229] DEBUG UFL_LEGACY: Blocks of each mode:
3 full
2023-10-28 16:40:29,650 [1229] DEBUG pulse.iterate: Solver did not converge...
2023-10-28 16:40:29,653 [1229] DEBUG pulse.iterate:
NOT CONVERGING
2023-10-28 16:40:29,653 [1229] DEBUG pulse.iterate: Reduce control step
2023-10-28 16:40:29,654 [1229] DEBUG pulse.iterate: Assign old state
2023-10-28 16:40:29,661 [1229] DEBUG pulse.iterate: Check target reached: NO
2023-10-28 16:40:29,662 [1229] DEBUG pulse.iterate: Maximum difference: 1.000e-01
2023-10-28 16:40:29,662 [1229] DEBUG pulse.iterate: Try new control
2023-10-28 16:40:29,663 [1229] DEBUG pulse.iterate: Current control: 0.010
*** Warning: PETSc SNES solver diverged in 2 iterations with divergence reason DIVERGED_FNORM_NAN.
2023-10-28 16:40:33,643 [1229] DEBUG pulse.iterate:
SUCCESFULL STEP:
2023-10-28 16:40:33,644 [1229] DEBUG pulse.iterate: Check target reached: NO
2023-10-28 16:40:33,645 [1229] DEBUG pulse.iterate: Maximum difference: 9.000e-02
2023-10-28 16:40:33,646 [1229] DEBUG pulse.iterate: Check target reached: NO
2023-10-28 16:40:33,647 [1229] DEBUG pulse.iterate: Maximum difference: 9.000e-02
2023-10-28 16:40:33,648 [1229] DEBUG pulse.iterate: Try new control
2023-10-28 16:40:33,649 [1229] DEBUG pulse.iterate: Current control: 0.020
2023-10-28 16:40:35,673 [1229] DEBUG pulse.iterate:
SUCCESFULL STEP:
2023-10-28 16:40:35,675 [1229] DEBUG pulse.iterate: Check target reached: NO
2023-10-28 16:40:35,676 [1229] DEBUG pulse.iterate: Maximum difference: 8.000e-02
2023-10-28 16:40:35,677 [1229] DEBUG pulse.iterate: Adapt step size. New step size: 0.015
2023-10-28 16:40:35,678 [1229] DEBUG pulse.iterate: Check target reached: NO
2023-10-28 16:40:35,679 [1229] DEBUG pulse.iterate: Maximum difference: 8.000e-02
2023-10-28 16:40:35,682 [1229] DEBUG pulse.iterate: Try new control
2023-10-28 16:40:35,682 [1229] DEBUG pulse.iterate: Current control: 0.035
2023-10-28 16:40:37,685 [1229] DEBUG pulse.iterate:
SUCCESFULL STEP:
2023-10-28 16:40:37,686 [1229] DEBUG pulse.iterate: Check target reached: NO
2023-10-28 16:40:37,688 [1229] DEBUG pulse.iterate: Maximum difference: 6.500e-02
2023-10-28 16:40:37,690 [1229] DEBUG pulse.iterate: Adapt step size. New step size: 0.022
2023-10-28 16:40:37,691 [1229] DEBUG pulse.iterate: Check target reached: NO
2023-10-28 16:40:37,692 [1229] DEBUG pulse.iterate: Maximum difference: 6.500e-02
2023-10-28 16:40:37,694 [1229] DEBUG pulse.iterate: Try new control
2023-10-28 16:40:37,695 [1229] DEBUG pulse.iterate: Current control: 0.058
2023-10-28 16:40:40,110 [1229] DEBUG pulse.iterate:
SUCCESFULL STEP:
2023-10-28 16:40:40,111 [1229] DEBUG pulse.iterate: Check target reached: NO
2023-10-28 16:40:40,112 [1229] DEBUG pulse.iterate: Maximum difference: 4.250e-02
2023-10-28 16:40:40,113 [1229] DEBUG pulse.iterate: Adapt step size. New step size: 0.034
2023-10-28 16:40:40,114 [1229] DEBUG pulse.iterate: Check target reached: NO
2023-10-28 16:40:40,115 [1229] DEBUG pulse.iterate: Maximum difference: 4.250e-02
2023-10-28 16:40:40,116 [1229] DEBUG pulse.iterate: Try new control
2023-10-28 16:40:40,117 [1229] DEBUG pulse.iterate: Current control: 0.091
2023-10-28 16:40:42,528 [1229] DEBUG pulse.iterate:
SUCCESFULL STEP:
2023-10-28 16:40:42,529 [1229] DEBUG pulse.iterate: Check target reached: NO
2023-10-28 16:40:42,530 [1229] DEBUG pulse.iterate: Maximum difference: 8.750e-03
2023-10-28 16:40:42,531 [1229] DEBUG pulse.iterate: Adapt step size. New step size: 0.051
2023-10-28 16:40:42,532 [1229] DEBUG pulse.iterate: Check target reached: NO
2023-10-28 16:40:42,533 [1229] DEBUG pulse.iterate: Maximum difference: 8.750e-03
2023-10-28 16:40:42,534 [1229] DEBUG pulse.iterate: Change step size for final iteration
2023-10-28 16:40:42,535 [1229] DEBUG pulse.iterate: Try new control
2023-10-28 16:40:42,535 [1229] DEBUG pulse.iterate: Current control: 0.100
2023-10-28 16:40:43,743 [1229] DEBUG pulse.iterate:
SUCCESFULL STEP:
2023-10-28 16:40:43,744 [1229] DEBUG pulse.iterate: Check target reached: YES!
2023-10-28 16:40:43,746 [1229] DEBUG pulse.iterate: Check target reached: YES!
([Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 47),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 175),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 224),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 273),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 330),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 387),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 420)],
[Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 45),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 173),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 222),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 271),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 328),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 385),
Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 418)])
# Get displacement and hydrostatic pressure
u, p = problem.state.split(deepcopy=True)
V = dolfin.VectorFunctionSpace(mesh, "CG", 1)
u_int = interpolate(u, V)
new_mesh = Mesh(mesh)
dolfin.ALE.move(new_mesh, u_int)
fig = plot(mesh, show=False)
fig.add_plot(plot(new_mesh, color="red", show=False))
fig.show()